Maps

In this notebook we will build maps for the species of Anchylorhynchus, using Natural Earth as base map and combining species with compatible distributions in a single map.

rm(list=ls())
library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──

## ✓ ggplot2 3.3.2     ✓ purrr   0.3.4
## ✓ tibble  3.0.1     ✓ dplyr   0.8.5
## ✓ tidyr   1.0.3     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(rnaturalearth)
library(raster)

## Loading required package: sp

## 
## Attaching package: 'raster'

## The following object is masked from 'package:dplyr':
## 
##     select

## The following object is masked from 'package:tidyr':
## 
##     extract

library(RStoolbox)
library(ggthemes)

## Warning: package 'ggthemes' was built under R version 4.0.1

library(sf)

## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0

library(gridExtra)

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:dplyr':
## 
##     combine

Let’s read specimen data, keep only georreferrenced data and remove duplicates:

sp_data = readxl::read_excel('especimes estudados_2020_rev2.xlsx')

sp_data = sp_data %>%
  filter(str_detect(updated_locality,'mistake',negate=T) | is.na(updated_locality)) %>%
  dplyr::select(accepted_name,latitude, longitude) %>%
  filter(!is.na(latitude)&!is.na(longitude)) %>%
  distinct %>%
  arrange(accepted_name, latitude, longitude) %>%
  rename(species = accepted_name)

sp_data

Now, let’s download the base map (commented because only run once).

#basemap = ne_download(scale = 'large', type = 'GRAY_HR_SR_OB_DR', category = 'raster', destdir = './map')

After map is downloaded, let’s crop it to the region of interest:

basemap = raster::stack('map/GRAY_HR_SR_OB_DR/GRAY_HR_SR_OB_DR.tif')
bbox = raster::extent(c(xmin=-95,xmax=-25,ymin=-45,ymax=19))
cropped = raster::crop(basemap,bbox)
map_layer = ggRGB(cropped, ggLayer = TRUE,r = 1,g=1,b=1)
#the following is just to visualize the layer
print(ggRGB(cropped, r = 1,g=1,b=1))

countries = ne_countries(scale=50, returnclass = 'sf') %>%
  st_transform(crs = 4326) %>%
  st_crop(c(xmin = bbox@xmin,
            xmax = bbox@xmax,
            ymin = bbox@ymin,
            ymax = bbox@ymax))

## although coordinates are longitude/latitude, st_intersection assumes that they are planar

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

Now, let’s load a table indicating which species should be plotted together.

sp_map = readr::read_csv('map_combination.csv')

## Parsed with column specification:
## cols(
##   species = col_character(),
##   map = col_double()
## )

sp_map

Now, let’s create a function that plots each map:

plot_map = function(spp){
  
  this_occ = sp_data %>% 
    filter(species %in% spp$species) %>%
    mutate(species = str_c('A. ',species)) %>%
    st_as_sf(coords=c('longitude', 'latitude'), crs=4326)
    
  
  ggplot(this_occ) +
    map_layer +
    geom_sf(color='black',fill=NA, size=0.2, data=countries) +
    geom_sf(aes(shape = species, fill = species), size = 1.5, show.legend = "point") +
    coord_sf(expand = F, )+
    scale_fill_brewer(type='qual') +
    scale_shape_manual(values = 21:24) +
    #guides(colour = guide_legend(override.aes = list(size=0.25))) +
    theme(
    axis.text = element_text(size=7), 
    axis.ticks = element_line(),
    #legend.text = element_text(face='italic',size=6,margin=margin(t=.1,r=0,b=.1,l=.1,unit='pt')),
    legend.text = element_text(face='italic',size=6),
    legend.margin = margin(t=0,r=2,b=2,l=2,unit='pt'),
    legend.justification = c(1,1),
    legend.position = c(1,1),
    legend.title= element_blank(),
    legend.box.margin = margin(t=0,r=0,b=0,l=0,unit='pt'),
    legend.background = element_rect(fill = 'white',colour = NA),
    legend.box.background = element_rect(fill = 'white',colour = NA),
    legend.key = element_rect(fill = 'white',colour = NA),
    legend.key.height = unit(8,'pt'),
    legend.key.width = unit(5,'pt')
    ) +
    xlim(bbox@xmin,bbox@xmax) +
    ylim(bbox@ymin,bbox@ymax) 
}

Now let’s plot all maps:

plots = split(sp_map,sp_map$map) %>%
  purrr::map(plot_map)

plots = 1:length(plots) %>%
  purrr::map(~plots[[.x]] + labs(tag=LETTERS[.x]))

plots 

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

And now arrange them in a page:

p = grid.arrange(arrangeGrob(grobs=plots, 
                         nrow = 3))

ggsave('maps_2.pdf', plot=p, device='pdf',width=168,height=200,units='mm',useDingbats=F)
ggsave('maps_2.eps', plot=p, device='eps',width=168,height=200,units='mm')

plot(p)